
### This script runs the analysis of the British Household Panel Data contained in the online appendix as reported in tables A15, A17 and A18.

### use Risk and self-respect dataverse download folder as your working directory. Where needed, we refer to this as "yourWorkingDirectory"

## set working directory
#setwd("yourWorkingDirectory")

### run merge and prep code
source("rsr_bhps_1loadprocess.R") ## as noted in the script, this file may take some time to run

## load libraries 
library(survival)
library(stargazer)
library(sandwich)
library(AER)

### linear probability model (OLS)
felmgour.3 = lm(rsrb2ss.demeaned ~ gour + sex +  marital.f + edqual + self.employed + age + age2 + log.income, data = long.ar[which(long.ar$risk.status1 == 0),])
rob.vc.felmgour.3 <- vcovHC(felmgour.3, type = "HC1")
ctfg.3 = coeftest(felmgour.3, rob.vc.felmgour.3)

### conditional logit model

long.ar$wave.fs = long.ar$wave.f
long.ar$wave.fs[which(long.ar$year < 2001)] = NA
long.ar$wave.fs = as.numeric(levels(long.ar$wave.fs))[long.ar$wave.fs]
long.ar$wave.fs = as.factor(long.ar$wave.fs)

cl.gour2.3e = clogit(rsrb2ss ~ gour + sex +  marital.f + edqual + self.employed + age + age2 + log.income + strata(pid), data = long.ar[which(long.ar$risk.status1 == 0),], method = "efron")

## creates coefficient estimates and model statistics for conditional logit model
stargazer(felmgour.3, cl.gour2.3e, style = "ajps", omit = c("sex", "marital", "edqual", "self", "age", "income", "Constant"), dep.var.labels = c("Demeaned binary SR", "Binary SR"), covariate.labels = "Gender-occupation unemployment risk", keep.stat = c("n", "rsq", "ll", "logrank"), digits = 3, digits.extra = 0, type = "text", title = "Table A15: Occupation-gender unemployment risk and self-respect in the BHPS, 2000-2008", notes.append = T, notes = c("Models include controls for gender, marital status, education, self-employment, income, age, and age squared.")
)


#### Table A17: The effect of overcoming adversity on self-respect. 
## Linear models of demeaned binary outcome, robust standard errors in parentheses. * indicates statistical significance at p < 0.05. All models include controls for educational achievement, income, marital status, self-employment, age, age2

all.cl = clogit(rsrb2ss ~ risk.status1 + bounded.prev.ar + marital.f + edqual.f + age + age2 + log.income + wave.f + strata(pid), method = "efron", data = long.ar)

all.lpm = lm(rsrb2ss.demeaned ~ bounded.prev.ar + risk.status1  + marital.f + edqual.f  + age + age2 + log.income + wave.f, data = long.ar)
rob.vc.lpm <- vcovHC(all.lpm, type = "HC1")
cft.lpm = coeftest(all.lpm, rob.vc.lpm)

## state specific demeaning
resolved.nar.rsrb2.1 = lm(rsrb2ss.demeaned.nar ~ bounded.prev.ar + risk.status1  + marital.f + edqual.f  + age + age2 + log.income + wave.f, data = long.ar)

resolved.ar.rsrb2.1 = lm(rsrb2ss.demeaned.ar ~ bounded.prev.ar + risk.status1  + marital.f + edqual.f  + age + age2 + log.income + wave.f, data = long.ar)

## state specific demeaning with interaction
resolved.nar.rsrb2.1i = lm(rsrb2ss.demeaned.nar ~ bounded.prev.ar*risk.status1  + marital.f + edqual.f  + age + age2 + log.income + wave.f, data = long.ar)

resolved.ar.rsrb2.1i = lm(rsrb2ss.demeaned.ar ~ bounded.prev.ar*risk.status1  + marital.f + edqual.f  + age + age2 + log.income + wave.f, data = long.ar)

## robust SEs and coefficient tests
rob.vc.nar.rsrb2.1 <- vcovHC(resolved.nar.rsrb2.1, type = "HC1")
cft.nar.rsrb2.1 = coeftest(resolved.nar.rsrb2.1, rob.vc.nar.rsrb2.1)
rob.vc.ar.rsrb2.1 <- vcovHC(resolved.ar.rsrb2.1, type = "HC1")
cft.ar.rsrb2.1 = coeftest(resolved.ar.rsrb2.1, rob.vc.ar.rsrb2.1)
rob.vc.nar.rsrb2.1i <- vcovHC(resolved.nar.rsrb2.1i, type = "HC1")
cft.nar.rsrb2.1i = coeftest(resolved.nar.rsrb2.1i, rob.vc.nar.rsrb2.1i)
rob.vc.ar.rsrb2.1i <- vcovHC(resolved.ar.rsrb2.1i, type = "HC1")
cft.ar.rsrb2.1i = coeftest(resolved.ar.rsrb2.1i, rob.vc.ar.rsrb2.1i)

stargazer(all.cl, all.lpm,  resolved.nar.rsrb2.1, resolved.ar.rsrb2.1, resolved.nar.rsrb2.1i, resolved.ar.rsrb2.1i,  omit = c("wave", "edqual", "marital", "age", "income", "Constant"), style = "ajps", keep.stat = c("rsq", "ll", "logrank"), dep.var.labels = c("Binary semi-strict", "", "","","","")
, covariate.labels = c("Current bad state", "Successful resolution","Successful resolution x current bad state"), column.separate = c(1,1, 1, 1, 1, 1), column.labels = c( "All", "All", "Good state", "Bad state", "Good state", "Bad state"), title = "Table A17: The effect of overcoming adversity on self-respect." , notes.append = T, notes = c("All models include controls for educational achievement, income, marital status, self-employment, age, age squared."), type = "text")


### Table A18: unemployment risk and overcoming adversity

felmgour.oa = lm(rsrb2ss.demeaned ~ gour + bounded.prev.ar + sex +  marital.f + edqual + self.employed + age + age2 + log.income, data = long.ar[which(long.ar$risk.status1 == 0),])
rob.vc.felmgour.oa <- vcovHC(felmgour.oa, type = "HC1")
ctfg.oa = coeftest(felmgour.oa, rob.vc.felmgour.oa)
ctfg.oa

cl.gour.oa = clogit(rsrb2ss ~ gour + bounded.prev.ar + sex +  marital.f + edqual + self.employed + age + age2 + log.income + strata(pid), data = long.ar[which(long.ar$risk.status1 == 0),], method = "efron")

stargazer(felmgour.oa, cl.gour.oa, style = "ajps", omit = c("sex", "marital", "edqual", "self", "age", "income", "Constant"), dep.var.labels = c("Demeaned binary SR", "Binary SR"), covariate.labels = c("Gender-occupation unemployment risk", "Successful resolution"), keep.stat = c("n", "rsq", "ll", "logrank"), digits = 3, digits.extra = 0, type = "text", title = "Table A18: Occupation-gender unemployment risk , overcoming adversity and self-respect in the BHPS, 2000-2008", notes.append = T, notes = c("Models include controls for gender, marital status, education, self-employment, income, age, and age squared.")
)
